## MODELING ENDOGENOUS EXPLANATORY VARIABLES

model {
	for (i in 1:N) {                     # N OBSERVATIONS
	    y[i] ~ dnorm (yhat[i], tau.y)
	    yhat[i] <- inprod(B[], X.1[i,]) + D[country[i],1]   
	    
	    X.1[i,5] ~ dnorm (x3.hat[i], tau.x3)
	    x3.hat[i] <- inprod(G[], X.2[i,]) + D[country[i],2]
	    
	    X.1[i,8] ~ dnorm (z3.hat[i], tau.z3)
	    z3.hat[i] <- inprod(H[], W.1[i,]) + D[country[i],3]	    
    }

	for (j in 1:J) {
		D[j,1:K] ~ dmnorm (mu.D[1:K], Tau.D[1:K,1:K])
	}

	for (k in 1:K) {
		mu.D[k] <- 0
	}
	
	# COVARIANCE MATRIX OF JOINT DISTRIBUTION
	Tau.D[1:K,1:K] ~ dwish (W[1:K,1:K], df)
#	df <- K+1
	
	Sigma.D[1:K,1:K] <- inverse(Tau.D[1:K,1:K])

	# COMPUTE THE VARIANCE TERMS	
	sigma2.d <- Sigma.D[1,1]
	sigma.d <- sqrt(sigma2.d)
	
	sigma2.x <- Sigma.D[2,2]+sigma.x3^2
	sigma.x <- sqrt(sigma2.x)

	sigma2.z <- Sigma.D[3,3]+sigma.z3^2
	sigma.z <- sqrt(sigma2.z)

    # DERIVE THE COVARIANCE	
	cov.dx <- Sigma.D[1,2]	
	cov.dz <- Sigma.D[1,3]
	
	# COMPUTE CORRELATION
	rho.dx <- cov.dx/sqrt(sigma2.d*sigma2.x)
	rho.dz <- cov.dz/sqrt(sigma2.d*sigma2.z)

    # PRIOR OF THE PRECISION PARAMETER
    tau.y ~ dgamma (0.01, 0.01)   # dgamma (r, mu), WHERE r IS THE SHAPE AND mu IS THE RATE PARAMETER      
    sigma.y <- 1/sqrt(tau.y) 
    
    tau.x3 ~ dgamma (0.01, 0.01)      
    sigma.x3 <- 1/sqrt(tau.x3)    
    
    tau.z3 ~ dgamma (0.01, 0.01)      
    sigma.z3 <- 1/sqrt(tau.z3)        
    
    # PRIORS ON UNMODELED COEFFICIENTS
    for (k in 1:M.1) {
    	B[k] ~ dnorm (0, 0.0001)
    }

    for (l in 1:L){
    	G[l] ~ dnorm (0, 0.0001)  
    }   

    for (q in 1:Q){
    	H[q] ~ dnorm (0, 0.0001)  
    }   
       
} # END OF MODEL


